A numerical study of a lifted H2/N2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {H}_2/\hbox {N}_2$$\end{document} flame excited by an axial and flapping forcing

The large eddy simulation method combined with the Eulerian stochastic field approach has been used to study excited lifted hydrogen flames in a stream of hot co-flow in a configuration closely corresponding to the so-called Cabra flame. The excitation is obtained by adding to an inlet velocity profile three types of forcing [(1) axial; (2) flapping; (3) combination of both] with the amplitude of 15% of the fuel jet velocity and the frequency corresponding to the Strouhal numbers St=0.30,0.45,0.60and0.75\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$St=0.30,\, 0.45,\, 0.60\, \,\text {and}\, \,0.75$$\end{document}. It is shown that such a type of forcing significantly changes the lift-off height (Lh\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_h$$\end{document}) of the flame and its global shape, resulting in the flame occupying a large volume or the flame, which transforms from the circular one into a quasi-planar one. Both the Lh\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_h$$\end{document} and size of the flames were found to be a function of the type of forcing and its frequency. The minimum value of Lh\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_h$$\end{document} has been found for the case when the axial and flapping forcing were combined and acted at the forcing frequency close to the preferred one in the non-excited configuration. The impact of the flapping forcing manifested through a widening of the flame in the flapping direction. It was shown that the excitation increases the level of temperature fluctuations caused by an intensified mixing process. The computational results are validated based on the solutions obtained for the non-excited flame for which experimental data are available.

www.nature.com/scientificreports/ stable combustion regimes and visibly shifts the blow-off limits. Research on an excited lifted non-premixed flame in a hysteresis regime, i.e., when depending on initial conditions a flame can be attached to a nozzle or remain lifted for the same fuel velocity, was performed by Demare and Baillot 16,17 . It was shown that by changing the amplitude and frequency of excitation one can enhance the combustion process or produce large fluctuations, and thus, weaken the flame stability. Kozlov et al. 18 analyzed micro-flames in the field of transverse sound waves. They found that the excitation can flatten the round flames and transform them into nearly plane flames. Surprisingly, for particular forcing frequencies the excitation led to a splitting of the flame into two separate branches in a very similar way as observed for bifurcating non-reacting, constant density jets 8 . More recently, the occurrence of the bifurcating phenomenon in flames was reported in numerical studies of Tyliszczak 19 focused on a low Reynolds number hydrogen flame. Application of only the flapping excitation caused the flame to change its initial circular shape into the planar one with two co-existing separate arms. Moreover, for some range of the excitation frequency, a triple-flame occurred. As discussed above, the excitation can alter the flame dynamics and influence pollution emission. In the present work we apply the LES method in the combination with the Eulerian Stochastic Field (ESF) approach 21 for the combustion modelling, and focus on the global impact of the excitation on the flame by applying three different excitation types (axial, flapping, axial plus flapping) with different forcing frequencies. Although the basic flow configuration closely corresponds to a well known Cabra flame 20 at the Reynolds number equal to 23600 the addition of the excitation is a novel element of the research. There are no experimental results for the excited Cabra flame, and hence, the credibility of the obtained results is proven by comparison with the measurements data available for the non-excited case. The present research is an exploratory numerical study in which we assess large-scale effects of the excitation, i.e., the change of the lift-off height or the change of the size and shape of the flame. It is shown for the first time that the dynamics of the lifted flame can be effectively controlled by a proper choice of the excitation type and its frequency.

Mathematical approach
We consider a low Mach number reacting flow for which the continuity, Navier-Stokes and transport equations of scalars within the framework of the LES method are defined as: where the bar and tilde symbols denote filtered quantities. The variables in Eqs. (1) are the velocity vector u , the density ρ , the hydrodynamic pressure p, the species mass fractions Y α and enthalpy h. The subscript α represents the index of the species α = 1, . . . , N-species . The quantities τ and D α , D are the viscous stress tensor and mass and heat diffusivities. The sub-filter tensor is given by τ SGS = 2µ t S , where S is the rate of strain tensor of the resolved velocity field and µ t is the sub-filter viscosity modelled as in 22 . The sub-filter diffusivities in the species and enthalpy transport equations are computed as D SGS = µ t /(ρσ ) where σ is the turbulent Schmidt or Prandtl number assumed equal to 0.7 23 . The set of equations (1) is complemented with the equation of state p 0 = ρR T with p 0 being the constant thermodynamic pressure and R is the gas constant.
The chemical source terms ρẇ α represent the net rate of formation and consumption of the species. A highly non-linear nature of this term means that sub-grid fluctuations cannot be ignored. In the present work, the scalar equations (species and enthalpy) are replaced by an equivalent evolution equation for the densityweighted filtered PDF function, which is solved using the stochastic field method proposed by Valiño 21 . Each scalar φ α is represented by 1 ≤ n ≤ N s stochastic fields ξ n α such that φ α = 1/N s N s n=1 ξ n α . The stochastic fields evolve according to: where the total diffusion coefficient is defined as Ŵ = D α + D SGS α , the micro-mixing time scale equals to τ SGS =ρ� 2 /(µ + µ t ) with � = Vol 1/3 cell being the LES filter width and dW represents a vector of Wiener process increments different for each field. Following Jones and Navarro 24 eight stochastic fields have been used. The test computations performed with sixteen fields did not show any substantial changes in the dynamics of the flame. Numerical method. We apply an in-house numerical code (SAILOR) based on the projection method for the pressure-velocity coupling 25 . The time integration is performed by means of an operator splitting approach where the transport in physical space and chemical terms are solved separately. The convective and diffusive parts of the governing equations are advanced in time using a predictor-corrector technique with the 2nd order Adams-Bashforth / Adams-Moulton methods. The chemical reactions are computed using a CHEMKIN interpreter. In the present study, we analyze the hydrogen combustion using a detailed mechanism of Mueller 26 involving 9 species and 21 reactions. The reaction terms are stiff and therefore they are integrated in the time applying the VODPK 27 solver that is well suited for stiff systems. The spatial discretization is performed on halfstaggered meshes applying the 6th order compact difference approximation for the Navier-Stokes and continuity equations 25,28 . The convective terms in the stochastic field equations are discretized applying TVD scheme with

Computational configuration
The basic flow configuration analyzed in this study corresponds to the so-called Cabra flame 20 , which we modify by adding the excitation at the nozzle exit. A schematic view of the computational geometry is shown in Fig. 1.
where u mean (x) is the mean velocity profile corresponding to the fully developed pipe flow (1/7 profile) and u turb (x, t) = 0.05U j represents turbulent fluctuations computed applying a digital filtering method proposed by Klein et al. 32 . This method guarantees properly correlated velocity fields which reflect realistic turbulent flow conditions. The forcing component u excit (x, t) is added to the streamwise velocity only and it is defined as: which is the superposition of axial and flapping forcing term with amplitudes A a and A f and frequencies f a and f f . The computations have been performed on three meshes consisting of 120 × 192 × 120 (coarse), 120 × 264 × 120 (medium) and 192 × 264 × 192 (dense) nodes compacted radially towards the flame region and axially towards the fuel nozzle. Preliminary tests showed that the medium mesh provides virtually the same solutions as the ones obtained on the dense mesh. This is due to the high-order discretisation method applied that leads to grid independent results already on relatively coarse meshes. However, to ensure a better resolution of the small scale phenomena the main computations have been performed on the dense mesh. Figure  As mentioned in the Introduction section, the simulations and results discussed in this study have an exploratory character as no experiments or numerical data are available for validation of the obtained results regarding the excited flames. Therefore the comparison was performed only based on measurements for the original Cabra flame configuration of which a spatiotemporal complexity is not much different from the cases with the excitation. As will be presented, the agreement between the present solutions and the experimental data is sufficiently good to assume that the obtained results are accurate and reliable. www.nature.com/scientificreports/

Results and discussion
Three-dimensional flame behavior. The fuel issuing from the nozzle mixes with the co-flowing hot stream and auto-ignites. The ignition spots appear at the locations of the most reactive mixture fraction ξ MR = 0.053 at distances far from the inlet, i.e. y/D ≈ 20 . Then, the flame spread radially, propagates upstream and stabilizes as a lifted flame in between y/D ≈ 7.5 − 11.0 depending on the test case. Figure 2 shows isosurfaces of an instantaneous temperature ( T = 1200 K ) and Q-parameter ( Q = 10 s −2 ) for the cases without the excitation and for A 45 , A 60 and A 75 . The Q-parameter is commonly used to indicate organized vortical motion. It is defined as Q = 1/2(� ij � ij − S ij S ij ) , where S ij and ij are the symmetric and antisymmetric parts of the velocity gradient tensor. Here, it visualizes the effect of the excitation, which manifests by the occurrence of toroidal vortices formed in the vicinity of the nozzle. They mutually interact through rib-like vortices and break up further downstream. Compared to the unexcited case (see Fig. 2a) the vortices are much more pronounced, the distances between them are dependent on the forcing frequency and decrease with increasing St a . One can notice that in the case A 60 the jet at y/D ≈ 1.0 − 4.0 is wider than for A 45 and A 75 and its shape reminds a barrel. The temperature distributions reveal that the flame lift-off height depends on the forcing frequency and turns out to be the smallest for A 60 that will be further confirmed by time-averaged results. The flame behaviour changes significantly when the flapping forcing is turned on. Figure 3 shows the results for the cases AF 60 and AF 75 in which the flapping forcing causes that the toroidal vortices are tilted about the axial direction. Note that this effect is noticeably stronger for AF 75 . For AF 60 the tilting of the vortical rings is overwhelmed by its very strong amplification by the axial forcing term, as will be discussed in the next sections. In both cases, however, the rings tend to move alternately to the left and right side of the domain and in effect, the flames with the flapping excitation become wider in the 'x'-direction and narrower in the 'z'-direction. In Tyliszczak A. 19 it was shown that for a low Reynolds number ( Re = 4000 ) the flame can even bifurcate (i.e. split into two separate branches), however, this phenomenon does not occur here. For the present case with Re = 23600 and a relatively high level of the inlet turbulence intensity the toroidal structures are destroyed before reaching a bifurcation point that usually is located at y/D ≈ 5 33 . www.nature.com/scientificreports/ Effect of the excitation in spectral space. Figure 4a shows an axial velocity spectrum computed based on a velocity signal recorded in the axis at y/D = 3.0 for the case without the excitation. It can be seen that there are no distinct peaks in the spectrum that could be related to the preferred mode frequency or the parring process. Apparently, as can be seen in Fig. 2a, the level of turbulence imposed at the inlet prevents the formation of strong, well defined vortical structures of which periodic occurrence would certainly manifest also in the spectrum. Instead, only a rise of the amplitudes of the fluctuations at a broadband range of the frequencies is observed, which is centred around St = 0.6 . The excitation at St a = 0.60 was chosen to match exactly this value, whereas the excitation at St a = 0.30 corresponds to its sub-harmonic at which the parring process could exist. Figure 4b shows the velocity spectra for the cases with the axial forcing only. Distinct peaks corresponding to the excitation frequencies are readily seen as they are pronouncedly larger than the background level. The cases A 45 and A 75 do not show anything exceptional. The peaks related to their basic frequencies are virtually the only ones visible. Further downstream they become wider and lower (not presented) and it seems that there are no additional phenomena created by the forcing at these frequencies. The results for A 30 and A 60 are significantly different. In the former case (green line in Fig. 4b) the excitation causes intensified velocity fluctuations not only at the basic frequency but also at its harmonic St = 0.60 . The most striking difference is, however, for the case A 60 for which the whole bench of highly energetic harmonics is found. They appear as the results of interactions between subsequent vortices. These interactions take place through the elongated rib structures (see Fig. 2c) that connect the vortical rings. One could expect that existence of harmonics causes an intensified mixing at small scales that speeds up the ignition process. The velocity spectra obtained for the cases with the axial and flapping forcing acting together show similar features, though, the harmonics are much weaker. The spectra for the cases with the only flapping excitation turned on are not significantly different from the 'no forcing' case as in the axis locations close to the inlet the impact of the flapping should not be pronounced by definition.

Impact of the excitation on the lift-off height.
The lift-off height ( L h ) of the flames is estimated based on the time-averaged results. The time-averaging procedure started when the flames were fully developed and it lasted for a period of at least 300D/U, which was found sufficient to obtain well convergent statistics. Figure 5 (a) AF 60 (b) AF 75  . It can be seen that both the shapes of the flames and L h are dependent on the type of the excitation. The L h was measured as the lowest point in the domain where the temperature or Y OH exceeded the above-defined thresholds. In all the cases these points occur not in the flame axis but a few diameters off-axis. The L h predicted based on the temperature criterion is slightly smaller than using the OH mass fraction criterion, however, the differences are not very significant ( �L h < 1.0D ). Worth noting is that both threshold lines predict very similar behaviour. In case A 60 in the central part of the flame these lines are almost straight and their inclinations to the flame axis depend on the forcing frequencies (not shown). When the flapping forcing is turned on the threshold lines become rounded. Figure 6 shows the PDF of the mixture fraction in the shear layer region (box of 5D × 1.5D × 1.5D located 2D below the flame front) for the cases with St a = 60 , and the dependence of the L h on the forcing frequencies for all analyzed configurations. A significant increase of the PDF around ξ ST = 0.476 for A 60 and AF 60 compared to the case F 60 and without the excitation means an intensified mixing process that should manifest by a shortening of L h . Indeed, for the case without the excitation L h = 9.7D ( L h = 10D in the experiment 20 ), which is visibly different from L h found when the excitation is applied. In these cases, depending on the forcing frequency and excitation type L h changes in the range 7.9D − 11.2D . Its lowest value occurs when the axial and flapping excitation act together. Turning on only the flapping mode causes that L h reaches the maximum at St a = 0.45 after which it continuously decreases as the effect of intensified mixing of smaller spatiotemporal flow scales caused by higher forcing frequencies. For the axial excitation, both acting solely or together with the flapping mode, L h behaves differently. It reaches the maximum at St a = 0.75 , but first, i.e., for St a < 0.75 , it suddenly drops down for St a = 0.45 and St a = 0.60 for the axial-flapping and axial modes, respectively. The occurrence of these minima can be to some extent related  www.nature.com/scientificreports/ the fuel burns in the shear layer ( y/D ≈ 11 ) and the near axis temperature distributions are correctly captured. Worth noticing is the fact that in the experiment the co-flow temperature was biased by the 3% error 20 that could lead to ±15 K difference compared to the assumed T cf = 1045 K . As observed by Navarro-Martinez and Kronenburg 35 the co-flow temperature has a substantial impact on the flame stabilization height. Nevertheless, it seems that ±15 K error in T cf has not a significant impact on the present results. Regarding the excited flames it can be readily seen that close to the inlet all profiles (also for the non-excitated case) are very similar and differences start to be seen only downstream. The temperatures are definitively larger in the axis for all excited cases, whereas for the cases AF their maxima move radially towards larger x/D locations. In general, the axial excitation causes a faster temperature rise, while the flapping forcing makes the temperature more uniform along the 'x'direction. The profiles along the 'z'-direction (Fig. 7b) show that further from the nozzle the flapping excitation can lead to 30% narrowing of the flame in respect to the 'x'-direction. Taking into account the fluctuations profiles one can see that the excitation rises their level close to the nozzle. Further downstream and in the vicinity of L h location the flapping mode flattens the profiles and shifts their maxima in the 'x'-direction. This is the effect of a more intense mixing with the co-flowing stream that makes the excited flames wider.

Conclusions
The paper presented the LES studies of the H 2 /N 2 flame excited with the axial and flapping forcing. The correctness of the results was confirmed by comparison with available experimental data for the unexcited case. It was found that the lift-off height of the flame, its size and shape can be altered in a wide range depending on the type of the excitation and its frequency. Compared to the unexcited flame the lift-off height can be increased or decreased. Its minimum value has been found for the case with the combination of both the axial and flapping forcing at the frequency close to the centre of the broadband frequency range regarded as the preferred one in the non-excited configuration. The impact of the flapping forcing manifested through a widening of the flame in the flapping direction. www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.